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SIMM1 


A  LINEARIZED  PARTICLE  CODE 


I.  INTRODUCTION 


Seif-pinched  relativistic  electron  beams  propagating  in  resistive  plasma 
are  subject  to  the  resistive  hose  instability,  essentially  a  growing  snakelike 
oscillation  which  is  often  observed  to  destroy  the  integrity  of  the  beam.^ 
For  this  instability,  the  linear  regime  (small  instability  amplitude)  is  of 
particular  interest,  because  the  mode  grows  convectively  backward  in  the  beam, 
and  therefore  it  is  possible  within  the  linear  regime  for  the  instability 
amplitude  at  any  given  location  in  the  beam  to  reach  a  maximum  level  and  then 
decay. 

A  variety  of  analytic  techniques  have  been  used  to  study  the  linearized 
hose  instability,  '  but  these  approaches  have  not  been  able  to  treat  an 

essential  feature  of  the  problem,  namely  that  self-pinched  beams  generally  are 
not  in  time-independent  equilibrium.  The  very  front  of  the  beam  head  is 

always  unpinched  and  therefore  continuously  erodes  away;  also  the  beam  radius 
increases  steadily  due  to  scattering  off  gas  molecules,  and  the  particle 

energy  decreases  and  spreads  due  to  various  loss  mechanisms.  Even  the  radial 
profile  of  the  beam  density  may  change  with  time  and  from  place  to  place  in 
the  beam. 

Two  approximate  numerical  models  have  also  been  developed  previously  to 

simulate  the  space-time  evolution  of  the  linearized  resistive  hose 

a  7 

instability.  *  In  these  treatments,  the  actual  particle  dynamics  is  replaced 
by  simplified  macroscopic  models,  which  are  carefully  chosen  to  preserve  a  key 
feature  of  the  dynamics,  the  spread  in  betatron  oscillation  frequency  among 
the  electrons  following  various  oscillation  orbits  in  the  pinched  beam.  (If 
this  feature  were  absent  from  the  model,  the  hose  instability  would 
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erroneously  become  absolute  rather  than  convective.]  The  macroscopic  models 


are  capable  of  treating  time-dependent  equilibria,  but  they  are  limited  to 


cases  where  the  equilibrium  beam  density  profile  n^(r,2,t)  has  the  self¬ 
similar  form  n^[r/a(z ,t)] ,  and  furthermore  the  accuracy  of  the  models  has  been 


demonstrated  only  for  cases  where  this  profile  has  the  Bennett  form,  n^  »  n^Q 
[l+r^/a^(z,t)]“^.  However,  the  actual  radial  density  profile  of  a  self- 


pinched  beam  may  be  hollow,  or  may  have  wings  of  significant  amplitude  at 


large  radius,  or  may  (on  the  other  hand)  be  cut-off  sharply  at  some  radius,  or 


may  vary  greatly  in  shape  from  place  to  place  in  the  beam.  These  types  of 


details  of  the  equilibrium,  which  cannot  be  represented  in  any  of  the  existing 


models,  can  strongly  affect  the  growth  of  the  hose  mode. 


The  technique  generally  used  to  study  instabilities  in  temporally  and 


spatially  varying  situations  is  particle  simulation,  but  in  practice  a 


standard  unlinearized  simulation  using  any  reasonable  number  of  simulation 


particles  (e.g.  106  or  even  more)  would  have  a  noise  level  large  enough  to 


mask,  or  misrepresent  the  instability  in  the  linear  regime.  To  overcome  these 


problems  and  to  provide  a  very  general  treatment  of  the  resistive  hose,  we 


have  developed  a  linearized  particle  simulation  code,  which  successfully 


determines  both  growth  rates  and  saturation  amplitudes  during  the  small- 


amplitude  stage  of  the  instability.  The  techniques  used,  which  build  upon 


the  earlier  work  of  Friedman,  Denavit  and  Sudan,10  may  also  be  applicable  to 


other  instabilities  which  grow  out  of  complicated  or  time-varying  equilibria. 


The  code,  SIMM1,  treats  the  zero-order  axisymmetric  beam  dynamics  (m  ■  0, 


where  the  azimuthal  dependence  is  Fourier  analyzed  into  e  dependence)  as  a 


group  of  axisymmetric  macroparticle  rings.  This  portion  of  the  code  is  named 


SIMMO  and  has  been  discussed  in  a  previous  publication.1^  The  m  ■  1  particle 


dynamics  associated  with  the  hose  perturbations  is  represented  by  linearized 


i  vi  vrM 


displacements  of  these  rings  from  the  beam  axis  and  by  m  ■  1  bunching  of  the 
charge  density  around  the  ring. 

A  large  number  of  these  rings  are  followed  in  time  as  they  interact  with 

each  other  and  with  the  background  gas.  The  beam  is  assumed  to  be  highly 

-2  2  2 

relativistic  so  that  quantities  of  the  order  of  y  =  1-v  /c  are  neglected. 
The  paraxial  approximation  is  also  made,  i.e.,  the  axial  velocities  of  the 
particles  are  taken  to  be  c,  and  the  perpendicular  velocities  are  assumed  to 
be  small.  The  dynamics  of  the  particles'  motions  are  determined  from 
solutions  of  Maxwell's  equations.  A  simplified  form  of  Maxwell's  equations  is 
appropriate  to  highly  relativistic,  paraxial  beams. **  The  equations  are 
presented  in  Appendix  A.  The  beam  travels  through  a  gas  which  is  assumed  to 
be  highly  collisional  and  so  dense  that  its  mass  motions  can  be  neglected  on 
the  beam  time  scale.  The  beam  partially  ionizes  the  gas  which,  in  turn, 
provides  charge  neutralization  and  allows  the  beam  to  self-pinch  in  its  own 
magnetic  field.  The  gas  Is  modeled  as  a  stationary  conducting  medium  with  a 
local  scalar  conductivity  which  is  determined  by  the  processes  of  beam-impact 
ionization,  avalanche  ionization  (when  there  are  large  beam-generated  electric 
fields),  and  recombination.  Conductivity  physics  can  be  modeled  at  various 
levels  of  sophistication.  The  model  used  in  SIMM1  is  discussed  in  Appendix 
A.  The  code,  then,  consists  of  axisymmetric  dynamics,  linearized  m  -  1, 
dynamics,  electromagnetic  field  generation,  and  interactions  with  the 
background  gas. 

The  linearization  procedure  presents  difficulties  which  must  be  dealt 
with  to  determine  the  long-time  behavior  of  the  system.  The  zero-order 
axisymmetric  motion  of  the  rings  occurs  under  the  influence  of  an  anharmonic 
central  potential,  so  the  period  of  their  motion  depends  on  the  oscillation 
amplitude.  Within  a  linearized  model,  this  effect  leads  to  secular  behavior 


of  the  m  ■  1  perturbations  to  the  orbits.  With  a  large  enough  number  of 
particles,  the  individual  orbit  secularities  would  average  out  to  yield  non¬ 
secular  macroscopic  quantities,  but  to  run  a  practical  simulation  some  method 
of  suppresssing  the  long-time  orbit  secularities  must  be  employed.  We  have 
developed  a  method  of  coarse  graining  to  deal  with  this  problem. 

We  will  discuss  the  linearization  process  in  section  II;  the  problem  of 
secular  behavior  and  its  resolution  by  coarse  graining  in  section  III. 
Details  are  given  in  Appendices  A  and  B. 


II.  CODE  STRUCTURE 

SIKM1  is  greatly  simplified  by  our  approximation  that  all  particles  have 
axial  velocity  c,  so  that  no  particle  overtaking  is  permitted  along  the  beam 
axis.  As  a  result,  we  can  consider  the  beam  to  be  made  up  of  a  series  of 
slices  of  particles  stacked  along  the  beam  axis  from  its  head  to  its  tail. 
The  simulation  can  be  thought  of  as  a  number  of  two-dimensional  (r,0) 
simulations,  with  information  carried  along  the  axis  by  the  fields  and 
conductivity. 

We  use  as  our  independent  variables  r,9,c  =  ct-z,  and  z  where  t  is  the 
time  and  ct  is  the  position  of  the  beam  head  in  space.  The  coordinate  Q  is  a 
measure  of  the  distance  along  the  axis  from  the  beam  head.  In  this  set  of 
coordinates,  z  acts  as  a  time-like  variable  with  regard  to  beam  dynamics, 
while  C  is  a  constant  of  the  motion  for  any  particle. 

A.  Particle  Dynamics 


The  beam  particle  dynamics  is  represented  in  terras  of  the  motion  of 
particle  rings  under  the  m  ■  0  and  m  *  1  forces.  Each  ring  represents  the  set 
of  all  particles  with  the  same  initial  value  of  r,c,  radial  velocity,  and 


angular  momentum.  It  Is  obvious  that  under  the  influence  of  an  axisymmetric 
force  the  m  *  0  motion  of  this  set  of  particles  is  such  that  it  remains  a 
ring,  with  radius  increasing  or  decreasing  in  time.  To  calculate  the  m  *  0 
dynamics  of  the  ring,  it  is  thus  sufficient  to  follow  the  dynamics  of  any  one 
particle  on  the  ring.  It  is  more  convenient  to  do  this  in  Cartesian 

i  i  > 

coordinates  in  the  i—  0  plane;  our  method  is  described  elsewhere. 

,  To  first  order  in  treating  the  m  *  1  perturbations,  each  particle  can  be 

assigned  a  small  time-dependent  displacement  from  its  equilibrium  orbit  in 
the  r-9  plane.  Friedman,  et  al.*®  have  shown  that  to  first  order  the  set  of 
particles  forming  any  given  ring  still  lie  on  a  circular  ring,  but  that  the 
center  of  the  ring  is  displaced  from  the  origin  by  a  small  vector  quantity, 
and  that  in  addition  the  charge  density  around  the  ring  is  no  longer  uniform, 

i(e-0  ) 

but  rather  is  (slightly)  bunched  as  e  .  The  linearized  m  *  1  dynamics 

of  all  the  particles  in  the  ring  is  thus  characterized  by  four  time-dependent 
quantities:  the  two  components  of  the  displacement  of  the  ring  center,  the 

bunching  axis  9q,  and  the  amplitude  of  the  bunching.  These  four  quantities 
can  be  represented  in  various  ways;  the  following  scheme  is  one  way  of  doing 

this  and  of  tying  the  dynamics  of  the  whole  ring  to  that  of  one  reference 

particle  on  it. 

Figure  1  shows  the  equilibrium  and  perturbed  positions  of  some  ring.  R 
and  P  are  two  particles  on  the  unperturbed  ring,  and  R"  and  P"  are  the 
perturbed  positions  of  these  particles.  We  represent  the  line  segment 
RR'  ■  Re(e(R)J.  The  quantity  e(R)  is  the  perturbation  of  the  ring  at  the 
point  R.  In  the  same  manner  the  line  segment  FP"  -  (e(P)J.  Note  that  the 

arc  RP  may  not  be  the  same  length  as  the  arc  R'P',  since  the  ring  of  charge  is 

bunched  as  well  as  displaced. 
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It  can  be  shown  that  if  we  choose  the  relationship  between  e(R)  and  e(P) 
to  be  given  by 

er(P)  -  er(R)eU, 

and 

e9(P)  -  eQ(R)eH, 

where  is  the  azimuthal  separation  of  P  and  R,  the  motions  of  the 
perturbations  are  consistent  with  the  dynamics  of  the  ra  -  1  fields  and  that  we 
only  need  to  follow  the  dynamics  of  the  perturbation  at  a  single  reference 
point  R.  The  two  complex  quantities  ef  and  thus  contain  all  the 
Information  about  the  bunching  and  displacement  of  the  ring. 

Although  we  have  chosen  a  complex  representation  of  the  perturbations, 
physical  quantities  must  be  real.  If  we  let 

-*■  -*■ 

♦  r  i 

e  -  er  +  i  e\ 

we  can  write 

PP-  -  Re(e(R)ei<^ J 

■  (e£  cos$  -  e*  sin$,  Eg  cos<fr  -  Eg  sin$). 

We  choose  to  identify  e*  -  e*(R),  e*1  -  -  e*(R),  e*  =  eJ(R)  and  Eg1  -  -  Eg(R). 

PP"  ■  (e*  cos$  +  e*1  sin$,  e*  cos*  +  e*1  sin$). 


Thus 


|U|Wt!C|^J<lJ^JLw>j 


yr,  r  *."  *-  *-t 


sa 


m 


The  particle  dynamics  are  computed  in  Cartesian  coordinates.  In  that 

coordinate  system,  we  need  to  follow  the  four  quantities  e*,  e*,  e**,  e**, 

x  y  x  y 

defined  by 


c1’11  -  e1’11  coset  +  t1'11  sina, 
r  x  y 


I, II  .  A  I, II  . 

-  e  sina  +  t  *  sina, 
x  y 


where  a  is  the  azimuthal  coordinate  of  the  reference  particle  R. 

The  electromagnetic  field  equations  for  the  conditions  of  an  ultra- 
relativistic  electron  beam  in  a  collisional  gas  were  derived  by  Lee^  and  are 
given  in  Appendix  A  along  with  a  simplified  conductivity  equation  which  models 
direct  ionization,  avalanche  ionization  and  recombination  of  the  background 


We  consider  only  beam  perturbations  linearly  polarized  along  the  y 
axis.  For  this  case  the  linearized  current  density  and  conductivity  are 


JL  -  J°  +  j}  sine, 
bz  bz  bz 


o  =  a0  +  sin6. 


The  only  nonzero  electric  and  magnetic  fields  are 


B°  +  Bg  sine, 


B*  cos6, 


>,<  /•  A  . *•  •  .V-iSVAVr.'  «>>  ■  ,  -  „* .  .Vv-.'  \*Vf' 


q  dz  Gx^  *  Er  sina  cos®  "  Eg  cosa  sinot 


“  c  B^  sina  cosa  -  c  Bg*  cosa  sina 


I  2  I  3E  3B° 

+  (e  cos  a  +  e  cosa  sina)  (•*— -  -  c  •*—--} 

x  V  v3r  at  / 


12  T 

+  ( e  sin  a  -  e 
x  v 


cosa  sina)  (E°/r  -  c  B°/r)  , 
7  r  rj 


ul  C 

o  d  ,  »II  ,  TI  2  12 

q  dz  6x  J  "  E-  cos  a  +  En  8in  c 


+  c  B*  sin2a  -  c  B*1  cos2a 


3E°  3B? 


II  2  TT 

+  ^ex  cos  a  +  ey  cosa  sinaj 


+  sin2a  -  e*1  cosa  sinaj  (E°/r  -  c  B°/rJ 


-  «"  sin'a  +  E*  cos2a 


+  c  B*  cos2a  -  c  B*1  sin2a 

r  9 


I  T  9  dt*  dBa 

+  Ccx  sina  cosa  +  ey  sin  aj  (•—  -  c  j^~) 


+  (e*  co82a  -  e*  cosa  sinaj  (E°/r  -  c  B^/rJ^ 


V.vv 


MU 


2 

ra  c  ,  TT 

-2-4-  (t  s11) 

q  dz  y 


E**  cosa  sina  -  E?  cosa  sina 
r  0 


-  c  sina  cosa  -  c  B^  cosa  sina 
r  o 


II  IT  2  3Er  3B0 

+  sina  cosa  +  e“  sin  aj  -  c 


+  (.e11  cos^a  -  e**  cosa  sinaj  (E°/r  -  c  B°/rJ, 


.  •  de 

where  e  =  -r—  . 

dz 

These  four  complicated  equations  form  the  basis  of  the  perturbation 
dynamics.  Although  the  total  beam  perturbation  is  along  the  y  axis,  the 
motion  of  any  single  ring  is  such  that  its  center  can  lie  at  any  point  off  the 
axis.  Only  in  adding  all  the  ring  positions  together  do  we  get  the  effects  of 
the  simplified  perturbation.  Otherwise  stated,  a  single  ring  has  angular 
momentum,  so  its  motion  does  not  in  general  have  the  same  symmetry  as  the 
driving  fields;  however  (it  is  tacitly  assumed)  the  distribution  over  angular 
momentum  is  even,  which  restores  the  symmetry  of  the  beam  as  a  whole. 

B.  Current  Density 

In  the  paraxial  approximation,  the  beam  current  density  component  J_ 
along  the  axis  is  proportional  to  the  charge  density  and  the  components  of  3 
perpendicular  to  the  axis  are  neglected.  The  charge  density  for  a  beam  made 
up  of  paint  particles  formed  into  rings  can  be  written  as 


Jz  <r'9>  “  ^ st*--rk.(*)js(B-ok(*>j/*, 

where  the  sum  is  over  the  rings  and  the  angle  $  is  the  azimuthal  angle  of  Fig. 
1.  The  represents  the  displacement  of  the  ring  from  the  axis  and 

9k($)  represents  charge  bunching. 

For  a  system  of  finite  sized  particles,  the  6-function  6(r-rjt(^)J  is 
replaced  by  a  particle  shape  function  S^r-r^C^)  J.  Within  the  framework  of 
the  linearization  procedure,  we  can  write 

rk(9)  ■  rk  +  £rk  eU> 

9k<«  *  9k  +  e9k  '“'V 

3S 

SrU-rk(«J  -  yr-r»)  -  *rk  e1*  -J  <r-r”), 

-  e9k  e^/r°  6'(0-o£). 


The  zeroth  and  first  order  current  densities  can  be  written 


v 


1  i9  qkc  o  e0k 

Jze  *  '  I  2w  Sr(r_rk)  "V  le 


,Vfr,  o.  1<9_ak) 

I  2nr  ,  o  <r  rk  rk  * 

9rk 


Finally,  we  can  write  the  current  density  as  a  real  function  consistent 


with  beam  displacements  in  the  y  direction  as 


J2  "  -  l  ifc  Sr (r"rk)  4”  I'  <•  €xk  sln“k  *  Eyk  “^““k 


*  <_  E*k  sl“\  +  Eyk  MS\)slnokJ 

q^c  3S  j  j. 

+  I  15?  7T  l(exk  cos\  +  Eyk  sl«ak)sinak 
3rk 

+  ^Exk  cos\  +  £yk  sino^cosaj. 


Since  Che  perturbation  current  density  contains  derivatives  of  the 
particle  shape,  a  shape  function  with  continuous  derivatives  is  necessary.  We 
have  chosen  a  quadratic  particle  interpolation  scheme  as  suggested  by 
Friedman,  et  al. 


III.  COARSE  GRAINING 

Although  formally  we  have  derived  a  set  of  equations  to  describe  the 
motions  of  the  perturbation  quantities,  the  nature  of  the  beam  particle  motion 
results  in  the  breakdown  of  the  linearization  procedure  under  some 
circumstances.  The  perturbation  of  each  individual  particle  ring  grows 
without  bound  even  In  the  absence  of  m  ■  I  forces.  To  visualize  this 
behavior,  we  consider  the  simpler  problem  of  motion  of  particles  in  a  fixed 
axisymmetric  potential.  (Since  information  only  propagates  backwards  in  a 
beam  moving  at  vz  *»  c,  this  is  equivalent  to  the  problem  of  the  motion  of  the 
first  perturbed  slice  of  the  beam,  which  is  subject  to  axisymmetric  fields 
determined  by  the  previous  unperturbed  beam  slices.) 

The  motion  of  a  beam  particle  in  the  plane  perpendicular  to  the  beam 
axis,  subject  to  an  anharmonic  potential  well,  is  quite  complicated.  In 
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general,  -the  particle  will  trace  out  an  unclosed  rosette  orbit.  The  essential 

feature,  however,  is  that  the  radial  oscillation  period  and  the  azimuthal 

drift  rate  depend  on  the  particle's  angular  momentum  and  the  amplitude  of  the 

orbit.  Because  of  this  dependence,  two  particles  whose  amplitudes  are 

initially  infinitesimally  different  will  separate  in  time  until  their 

separation  becomes  comparable  to  the  particles'  amplitudes  themselves.  If  we 

think  of  the  difference  in  the  amplitudes  as  the  perturbation  quantity,  we  see 

that  there  is  secular  growth  of  the  perturbation,  which  is  proportional  to 

time  as  long  as  the  perturbation  remains  small  compared  to  the  orbit 

amplitude.  We  expect  that  for  a  beam  slice  subject  only  to  an  axlsymmetrlc 

force,  but  initially  displaced  from  the  axis,  phase  mixing  among  particles  of 

different  oscillation  frequency  will  eventually  return  the  slice  to 

axisymmetry.  In  an  idealized  simulation  with  an  infinite  number  of  particles, 

this  process  could  be  well  modeled  for  an  indefinite  time;  the  growing 

individual  particle  perturbations  would  average  out  when  macroscopic 

quantities  were  calculated.  However,  in  an  actual  simulation  with  N  particles 

per  slice,  the  noise  level  (for  quantities  like  the  mean  displacement  of  a 

—  1/2  - 

slice)  is  of  the  order  of  eN  ,  where  e  is  a  typical  perturbation  amplitude 
for  an  individual  particle.  Since  e  «  t,  the  noise  eventually  comes  to 
dominate  at  a  time  t  which  scales  as  Thus  Increasing  N  delays  the 

breakdown  of  the  linearized  model,  but  eventually  the  phase  mixing  process  is 


overwhelmed  by  a  growing  noise  signal.  We  illustrate  this  behavior  by  running 
a  simulation  of  a  single  slice  in  which  the  magnetic  pinch  force  Is  calculated 


from  Amperes  law  for  a  fixed  Bennett  distribution  of  current,  i.e., 

-2 

J  (r)  *  J  (I+r2/a2)  ,  giving  a  force  F(r)  -  -  e2J„r(l+r2/a2)-1.  The  beam 


o  o  '  °  ^  o 

profile  initially  is  taken  to  be  a  similar  Bennett  profile  displaced  from  the 

axis.  Figure  2  shows  the  time  behavior  of  the  beam  centroid  y,  for  two  cases 
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which  differ  only  in  the  number  N  of  simulation  particles.  As  N  increases, 
the  time  that  the  simulation  matches  the  correct  physics  (steadily  decaying 
oscillations  of  y)  increases,  but  in  all  cases,  the  noise  dominates  the 
simulation  at  late  enough  times. 

Not  all  particle  orbits  exhibit  secular  behavior  at  the  same  rate.  A 
Bennett-like  potential  is  nearly  harmonic  close  to  the  beam  axis.  Particles 
whose  orbits  lie  in  this  region  show  very  slow  secular  behavior  while 
particles  which  make  large  excursions  in  radius,  i.e.  those  with  high 
perpendicular  energy  but  small  angular  momentum,  show  the  fastest  secular 
growth  in  a  state  potential. 

For  the  hose  instability,  which  is  convective,  the  macroscopic  beam 
perturbations  at  any  given  location  C  in  the  beam;  grow,  saturate  and  then 
decay.  The  growth  phase  is  well  modeled  in  a  linearized  simulation  of  the 
type  we  are  discussing,  since  the  exponential  growth  of  the  ordered  motion 
easily  overwhelms  the  secular  growth  of  the  noise.  However  the  saturation  and 
decay  cannot  be  modeled  by  our  straightforward  linearized  simulation  unless  an 
extremely  large  number  of  particles  are  employed.  If  the  simulation  is  to  be 
useful,  some  method  of  eliminating  the  effect  of  the  secularitles  must  be 
used. 

In  order  to  eliminate  the  secular  growth  of  the  noise,  it  suffices  to 

break  up  the  ordered  secular  growth  of  the  individual  particle 

perturbations.  In  general,  we  have  no  Interest  in  individual  particle 

dynamics,  but  only  in  the  various  macroscopic  moments.  Furthermore,  every 

particle  ring  carries  with  it  both  axisymmetric  attributes  —  x,y,v  ,v  ,y  — 

x  y 

that  determine  the  position  and  motion  of  the  ring  and  the  transport  of 

information  from  one  radial  region  to  another  by  means  of  large  scale  particle 

,  .  I  II  I  II  »I  »II  *1  *11 

orbits,  and  linearized  attributes  —  e  ,  e  ,e,e  ,e,e  ,e,e  — 

’  x’  x  ’  y’  y  ’  x’  x  y  y 


Chat  carry  the  perturbation  information.  The  problem  of  secularities  affects 

only  the  linearized  quantities.  Our  task  is  thus  to  find  some  technique  that 

preserves  the  unperturbed  orbits  and  the  present  value  of  the  macroscopic 

quantities,  but  somehow  averages  out  the  perturbed  quantities  pertaining  to 

the  individual  particles  before  they  have  the  opportunity  to  grow  intolerably 

large.  In  particular,  the  beam  property  that  acts  as  the  source  for  the 

electromagnetic  field  equations  is  the  axial  component  J  (r,?,z)  of  current 

z 

density,  which  is  equal  to  cp(r,?,z),  where  p  is  the  charge  density,  since  all 
particles  have  vz  =  c. 

We  have  developed  a  technique  that  we  call  "coarse  graining"  that 
accomplishes  these  objectives.  Basically,  this  consists  of  periodically 
performing  averages  of  all  the  perturbed  quantities  e  and  e  over  all  the 
particles  in  a  given  cell  in  configuration  space.  The  e's  and  e's  for  all  of 
these  particles  are  then  reset  to  the  average  value  at  that  cell.  The 
averages  are  weighted  in  such  a  way  that  the  perturbed  current  density 
J^Cr.C.z)  is  left  unchanged.  The  unperturbed  orbits  are  not  altered  in  any 
way.  The  process  is  similar  in  spirit  to  the  coarse  graining  technique  widely 
used  in  numerical  solutions  of  the  Vlasov  equation,^  in  that  numerically 
troublesome  fine  structure  is  periodically  smoothed  out,  by  periodically  re¬ 
assigning  locally- smoothed  values  to  the  distribution  function  while 
preserving  physically  meaningful  large-scale  quantities.  The  algorithm  for 
performing  the  coarse  graining  is  presented  in  detail  in  Appendix  B. 

The  coarse  graining  should  be  performed  at  intervals  which  are  long 
compared  to  the  individual  simulation  time  steps,  but  short  compared  to  the 
times  for  secular  growth.  Since  the  oscillation  period  of  typical  particles 
in  the  axisymmetric  potential  well  (the  "betatron  period")  is  the 
characteristic  time  for  doubling  of  the  initial  individual  particle  e's,  it  is 


suggested  that  the  coarse  graining  interval  be  of  the  order  of  the  betatron 
period.  Of  course  the  macroscopic  time  evolution  must  be  insensitive  to  the 
exact  choice  of  coarse  graining  interval,  if  the  results  are  to  be  meaningful. 

To  illustrate  the  technique,  the  problem  of  a  single  displaced  beam  slice 
subject  to  a  fixed  Bennett  potential  —  Fig.  2  —  was  rerun  on  the  code  with 
coarse  graining.  The  averaging  process  was  carried  out  at  various 
intervals.  For  intervals  between  about  0.5  to  1.5  betatron  periods,  the 
results  were  insensitive  to  the  averaging  frequency,  as  shown  in  Fig.  3.  The 
displacement  dies  out  in  a  manner  similar  to  that  without  coarse  graining,  but 
the  late-time  noise  which  was  seen  in  the  previous  figure  does  not  reappear. 

For  static  potentials,  the  coarse  graining  algorithm  we  have  described 
works  quite  well.  When  the  potential  varies  in  time,  as  it  does  in  the 
problems  of  interest,  where  beam  perturbations  self-consistently  drive  the 
fields  for  instability  modes  like  the  resistive  hose,  another  effect  can  prove 
troublesome  numerically.  Chambers^  and  Friedman,  et  al.^  have  shown  that  a 
few  particles  experience  almost  exponential  growth  of  their  perturbed 
quantities  e,  in  a  stochastic  way  that  is  nearly  impossible  to  resolve  in 
simulations.  Furthermore,  this  effect  can  be  spread  to  a  large  number  of 
particles  by  the  coarse  graining  algorithm.  We  deal  with  this  problem  in  an 
ad  hoc  way  which  seems  to  have  only  a  minimal  effect  on  the  system  as  a 
whole.  At  each  time  step  we  determine  which  ring  in  the  slice  has  the  largest 
value  of  e.  We  replace  the  e  and  e  of  this  ring  with  values  averaged  over  all 
the  particles  in  the  slice.  This  replacement  tends  to  make  the  slice  behave 
slightly  more  like  a  rigid  beam,  but  at  each  time  step  affects  only  one 
particle  in  about  2000. 

The  effect  of  applying  both  algorithms  is  a  linearized  particle  code 
which  can  treat  the  hose  instability  in  a  way  which  accurately  reproduces  the 


Fig.  3  Mean  displacement  of  a  beam  slice  in  a  fixed  Bennett  potential 
as  a  function  of  propagation  distance.  The  coarse  graining 
algorithm  was  employed.  Each  slice  contained  1000  rings. 

Fig.  3a:  The  coarse  graining  was  applied  every  10  time  steps. 

Fig.  3b:  The  coarse  graining  was  applied  every  20  time  steps. 

Fig.  3c:  The  coarse  graining  was  applied  every  30  time  steps. 
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growth,  saturation  and  decay.  Figure  4  shows  the  results  of  a  hose  simulation 
with  SIMM1,  including  coarse  graining.  We  note  that  at  any  given  location 
Z  in  the  beam,  the  mode  grows  (in  good  agreement  with  linear  theory^) 
saturates  and  then  decays;  as  a  function  of  z,  the  mode  continues  to  grow 
convectively  backward  in  ?.  Figure  5  shows  the  displacement  of  a  single  slice 
of  the  case  of  Fig.  4.  Figure  6  shows  the  results  of  a  similar  run  without 
coarse  graining.  For  this  case,  the  mode  grows  as  before  but  after  it 

saturates,  the  secular  behavior  becomes  dominant  and  the  displacements 
continue  to  grow  without  bound. 
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Fig.  6 


Hose  displacement,  y,  as  a  function  of  propagation  distance,  z, 
for  a  beam  slice  of  a  simulation  similar  to  that  of  Fig.  4,  but 
without  coarse  graining  algorithm.  The  slice  was  located  50  cm 
from  the  beam  head.  The  hose  amplitude  grows,  saturates,  and 
then  regrows  due  to  secular  growth  effects. 
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Appendix  A 


The  field  equations  used  in  SIMM1  are  the  reduced  set  of  Maxwell's 
equations  derived  by  Lee.^  Using  the  paraxial  approximation  and  Doppler 
shifted  variables,  we  can  write 


V2  A 
1  z 


-  — j  , 

c  z  ' 


^  “  4irp  , 


WO  * 

n  t  Z  ^  3cp 

vi  *  'i  "  Jz~  W  "  0  * 


7  «  A  -  .  Z  +  i-X  at  0 

'l  '1  TT  3c  ’ 


where  A  is  the  vector  potential,  <p  is  the  scaler  potential,  the  subscript  z  is 
associated  with  the  axial  direction  and  the  subscript  1  is  associated  with  the 
radial  variables  r  and  6. 

The  electromagnetic  fields  can  be  obtained  from 
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J  “  J,  +  a  E  . 
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The  field  equations  simplify  further  if  we  drop  the  transverse  displacement 
current  and  write 
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Then  the  field  equations  can  be  written  as 


and 
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where 


Cl  -  A  -  4>. 
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We  introduce  cylindrical  coordinates,  and  make  a  linearized  multiple 
decomposition  of  the  aximuthal  dependence.  Then 


J.  -  J°  +  J?  sin6, 
bz  b  bo  * 


IflitJWASJ  ^.*^.'1  LUJ.UAM  MP1IJ.IJI  wm  IJJJJJ  I- 


CL-  CL°  +  ol  sin8, 


<|>  -  $°  +  4*1  sind, 


if  we  consider  only  displacements  along  the  y  axis. 
The  field  equations  are 
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The  code  determines  the  conductivity  from  the  linearized  version  of 


If  "  +  Yi  c  “  0r  °2’ 


where,  if  the  propagation  medium  is  air,  we  use  the  rough  model 
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<  -  8.5x10“*  cm/statcoulomb. 


.kl JWl 


i  a-*.  »g  tj  fc*  U'l 


where 


8f  ■  7.0x10  sec/ cm, 


y±  -  [ApS3/(1  +  BS  +  CS2  +  DS3) Jsec”1 , 


pi 

5 


S  =  E  2/p2, 


E  is  the  electric  field  in  statvolts/cm,  p  is  the  air  density  in  atmospheres. 


A  -  1.42xl0~4 
B  -  9.18xl0-6 


C 

i 


C  »  2.66xl0~lu 
D  -  2.82xl0“17. 


One  could  envision  improving  this  model  in  a  variety  of  ways,  such  as  making 
k  and  &r  temperature-dependent,  including  attachment,  including  the  effects  of 
various  minority  species  and  chemical  changes  in  complex  gases  such  as  air, 
and  modeling  nonlocal  or  magnetic  field  dependent  conductivity  if  the  electron 
mean  free  path  is  long  compared  to  macroscopic  lengths  or  Larmor  radii. 
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Appendix  B 


Consider  the  equation  for  (r,?)  rewritten  in  the  following  way 


JzI(r‘°  "  7^  ^  ^  Sr(rk;  r_rk) 

2  2irr  N  kfN  r  *  K 

s  s 
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where  the  first  summation  is  over  all  particles  in  cell  Ng  and  the  second 
summation  is  over  cells. 

The  two  quantities  which  must  be  preserved  by  the  averaging  process  are 
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Recall  that,  in  the  average,  only  the  e's  are  changed.  The  quantities  and 

S  (rP,  r-r?)  depend  on  axisymmetric  variables  and  are  the  same  before  and 
r  k  tc 

after  the  averages.  If  Fkl  and  F^2  are  preserved  in  the  average  so  are  linear 
combinations  of  them.  In  particular 


FU  +  Fkl  '  ^  +  *£>  sin2\  *  Ux  -  *?  cos2ok 


will  be  unchanged.  If  we  choose  to  average  in  such  a  way  that  ej£  +  , 

(£  -  e£  )  cosZa^,  sin2ak>  and  sin2ak  are  preserved  then  F^j  and  F 

will  be  unchanged. 
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then  the  averages  can  be  written  as 
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Finally,  we  associate  the  averages  with  the  e's  of  individual  rings  by 


using  the  particle  shape  functions. 

In  addition  to  averaging  over  the  e's,  some  method  must  be  employed  to 

average  the  e's.  The  method  which  suggests  itself  is  to  form  averages  which 

preserve  J These  averages  depend  not  only  on  the  particle  shapes  and 

positions  but  also  on  the  rotation  velocity  of  the  rings.  In  practice,  these 

averages,  especially  at  small  radii,  tend  to  be  dominated  by  a  few  rapidly 

spinning  rings.  As  a  compromise,  we  have  chosen  to  perform  unweighted 
•I  »I  *11  »II 

averages  of  e  ,  e  ,  e  ,  and  e  •  Me  have  performed  a  number  of  tests  on  the 
x  y  x  y 

e  averages  and  the  method  appears  to  be  satisfactory. 
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